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Abstract 

The Hasegawa-Wakatani equations, coupling plasma density and electrostatic potential through 
an approximation to the physics of parallel electron motions, are a simple model that describes 
resistive drift wave turbulence. We present numerical analyses of bifurcation phenomena in the 
model that provide new insights into the interactions between turbulence and zonal flows in the 
tokamak plasma edge region. The simulation results show a regime where, after an initial transient, 
drift wave turbulence is suppressed through zonal flow generation. As a parameter controlling the 
strength of the turbulence is tuned, this zonal flow dominated state is rapidly destroyed and 
a turbulence-dominated state re-emerges. The transition is explained in terms of the Kelvin- 
Helmholtz stability of zonal flows. This is the first observation of an upshift of turbulence onset 
in the resistive drift wave system, which is analogous to the well-known Dimits shift in turbulence 
driven by ion temperature gradients. 
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I. INTRODUCTION 



Fusion plasmas and other turbulent flows in quasi-two-dimensional (2D) geometry can 
undergo spontaneous transitions to a turbulence- suppressed regime. In plasmas they are 
known as L-H (low-to-high confinement) transitions and are studied intensively because 
they effectively enhance the confinement, through suppression of anomalous or turbulent 
particle and heat fluxes. It is now widely accepted that emergent zonal flows are crucial 
to achieving confinement improvement The L-H transition is associated with nonlin- 
early self-generated poloidal E x B shear or zonal flows |2J in the tokamak edge region, 
which comprises the transition zone from inner hot core plasma to the outer cold scrape-off 
layer. Zonal flows reduce anomalous transport by absorbing energy from drift waves and 
by shearing apart eddies which mediate turbulent transport, and thus play a key role in its 
regulation. 

In this paper we present the results of analytic and numerical investigations of transitions 
between turbulence-dominated and zonal-flow-dominated regimes, using the Hasegawa- 
Wakatani (HW) model jjj, 4] for electrostatic resistive drift wave turbulence in 2D slab 
geometry. We find that bifurcations in the model correspond to the onset of drift wave 
turbulence, the generation of zonal flows, and the re-emergence of turbulence as the zonal 
flows become unstable, and observe that this is drift wave turbulence analog of the Dimits 
shift [5j in ion temperature gradient turbulence. 

Three energetic subsystems interact to produce the complexity observed in L-H tran- 
sition dynamics: the kinetic energy of turbulence, the kinetic energy of shear flows, and 
the potential energy contained in density or pressure gradients. The three major governing 
processes are generation of turbulence by drift waves, self-organization of zonal flows, and 
destabilization of the zonal flows. The instabilities that lead to these changes correspond 
to bifurcations of equilibrium solutions of model equations. If a tunable parameter crosses 
a stability threshold the qualitative nature of the solution changes. We say that a primary 
instability occurs at a linear stability threshold of the equilibrium with zero background 



flow, which physically corresponds to the onset and growth of drift waves. Theoretical 
and experimental [7] studies have indicated that the generation of drift wave turbulence in 
plasmas may occur by the Ruelle-Takens mechanism [8|, in which a limit cycle generated by 
a Hopf bifurcation undergoes a Niemark-Sacker bifurcation to a torus, which may undergo 
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FIG. 1: Primary instabilities generate turbulence from a potential energy reservoir, secondary 
instabilities lead to the growth of shear or zonal flows at the expense of turbulence kinetic energy, 
and tertiary instabilities may destabilize the shear or zonal flows. Zigzag green arrows represent 
dissipative channels. 

one or more bifurcations to higher- dimensional tori before the motion becomes chaotic. 

However, to complicate this generic turbulence onset scenario, in plasmas zonal flows 
will be generated beyond the primary threshold due to an instability of the drift waves, 
effectively suppressing drift wave activity. This instability causing the zonal flow onset is 
termed a secondary instability. We can consider the turbulence to be well-developed at the 
secondary instability; i.e., for heuristic purposes we assume the Ruelle-Takens sequence to 
have already occurred. 

A strong candidate for this secondary instability mechanism is modulational instability 
{9I, [h]], a special case of nonlinear mode coupling whereby modulation of a small scale 
monochromatic wave can transfer energy non-locally to a longer wavelength structure due 
to the ponderomotive force effect leading to excitation of zonal flows. One might also expect 
an inverse energy cascade, endemic to quasi two-dimensional flows in general, whereby local 
mode coupling channels energy into large scale structures. 

A different mechanism for this secondary instability that generates zonal flows is Kelvin- 



Helmholtz (KH) instability [111 . |l2j]. In this scenario the KH instability may be driven 
by radially elongated drift wave eigenmodes. The KH mode of the drift waves necessarily 
possess a zonal flow component, and provide a natural mechanism for the zonal flow growth. 

As the zonal flows become more energetic they are subject to tertiary instability which 
breaks up the coherent zonal structuring of the flow into turbulent small scale eddies via KH 
instabilities of the zonal flows. The small scale turbulence may again cohere via secondary 
instabilities. These interactions are schematized in Fig. HJ 
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Nonlinear interactions between zonal flows and drift waves results in an upshift of the 
boundary in parameter space for the tertiary onset of turbulence. This is known as the Dim- 
its shift in ion-temperature-gradient (ITG) driven turbulence, and the turbulence suppressed 
regime was mapped by gyrokinetic and gyrofluid simulations {jj. 

The simplest approach that captures the essential physics underlying the problem is low- 
dimensional dynamical modeling and analysis 13|, [14], |l5j, [l6j, which can provide a very 



economical tool to predict the transition. However, the tradeoff with such highly coarse- 
grained modeling is that it necessarily whites out information, and may therefore miss 
important physics and predict unphysical singular behavior 15] . Thus we require validation 
of the low-dimensional modeling results by computational simulations of finer models. 
The HW model jj], 4] was developed to investigate anomalous edge transport due to 
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211 ] . It includes the 



collisional drift waves, and has been widely studied 
effects of inhomogeneous background density and parallel electron dynamics described by 
Ohm's law. The density gradient drives the drift waves, which are destabilized by the par- 
allel electron resistivity. Convective nonlinearity regulates the linear growth of the resistive 
drift wave instability, and a quasi-stationary state is achieved where the resistive coupling 
balances the input. The HW model is particularly simple yet includes the essential physics 
for studying the self-consistent generation of turbulence and growth and decay of coherent 



macroscopic structures such as zonal flows 



171 ]. even though it does not describe physics 



that can be important in specific situations, such as magnetic curvature, magnetic shear, 
and electromagnetic effects. 

We emphasize that the parallel electron motion is important for generation, stabilization, 
and destabilization of zonal flows. The parallel electron response given by the generalized 
Ohm's law leads to resistive coupling between the electrostatic potential and the de nsity 
fluctuations. In toroidal geometry this coupling does not act on the flux-averaged parts [22] , 
and in the original or unmodified HW model we do not observe zonal flows. Modification 
of the resistive coupling term, described in Sec. HIl enables the generation of zonal flows. 
This corresponds to the difference between the ITG and the ETG (electron-temperature- 
gradient) cases discussed by Jenko et al. 12j, who found that suppression of the secondary 
KH instability in the ETG case, due to the adiabatic electron response, is removed in the 
ITG limit. 

In Sec. HIl we describe the HW model and discuss the treatment of parallel electron 
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motions. Linear stability analysis of the zero-flow background is also given to calculate 
transition points in parameter space. Numerical simulation results are given in Sec. IIHI We 
carry out a systematic parameter survey to locate the transition from a zonal-flow-dominated 
state to a turbulent state. To examine the hypothesis that this transition may be ascribed 
to the tertiary KH instability of the zonal flow, we study the KH stability of the generated 
zonal flows in the HW model in Sec. IIVI and compare the KH stability threshold with the 
transition boundary determined by simulation. Discussions and conclusions are presented 
in Sec. [V] 



II. MODIFIED HASEGAWA WAKATANI MODEL 



The physical setting of the HW model may be considered as the edge region of a tokamak 
plasma of nonuniform density n = no(x) and in a constant equilibrium magnetic field B = 
BqVz. Following the drift wave ordering 23|, the ion vorticity ( = V 2 <p (<p is the electrostatic 
potential, V 2 = d 2 /dx 2 + d 2 /dy 2 is the 2D Laplacian) and the density fluctuations n are 
governed by the equations 

^C + MC} = «(^-n)- J DV 4 C, (i) 

^n + {<p,n} = a(<p-n)-K^--D\7\ (2) 

where {a, b} = (da/dx)(db/dy) — (da/dy)(db/dx) is the Poisson bracket, D is the dissi- 
pation coefficient. The background density is assumed to have an unchanging exponential 
profile: k = (d/dx) In no- Electron parallel motion is determined by Ohm's law with electron 
pressure p e = nT e , 

j z = -env ez = — [<p Inn , (3) 

rjdz \ e J 

assuming electron temperature T e to be constant (isothermal electron fluid). This re- 
lation gives the coupling between ( and n through the adiabaticity operator a = 
—T e /(rjn u) C ie 2 )d 2 /dz 2 appearing in Eqs. ([T|) and ([2]). In our 2D setting a becomes a con- 
stant coefficient when acting on the drift wave components of if and n by the replacement 
d/dz — > ik z , where 2n/k z = L\\ ^> L y is a length characteristic of the drift waves' phase 
variation along the field lines. However, for the zonal flow components, this resistive cou- 
pling term must be treated carefully because zonal components of fluctuations (k y = k z = 
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modes) do not contribute to the parallel current [22|]. Recalling that turbulence in the toka- 
mak edge region, where there is strong magnetic shear, is considered here, k y = should 
always coincide with k z = because any potential fluctuation on the flux surface is neutral- 
ized by parallel electron motion. Let us define zonal and non-zonal components of a variable 

/ as 

zonal: (/) = -!- J fdy, non-zonal: / = /-(/), 

where L y is the periodic length in y, and remove the contribution by the zonal components 
in the resistive coupling term in Eqs. (J]]) and (j2J). Subtraction of the zonal components from 
the resistive coupling term a(<p — n) — > a{(p — n) yields the modified HW (MHW) equations, 

^-( + { ip X} = a(<p-n)-DV 4 ( } (4) 
d dip 

—n + {ip,n} = a((p-n)-K—-DV 4 n. (5) 

Evolutions of the zonal components can be extracted from Eqs. (jlj) and (jSJ) by averaging in 
the y direction: 

|</> + | = 

where / stands for ( and n. 

Wakatani and Hasegawa found 4j that excitations of waves having k z that maximizes the 
linear growth rate (for given k x and k y ) are most likely to occur, since the plasma can choose 
any parallel wavenumber (k z ). Using the parallel wave number of the maximum growth rate, 
a is given by a = Ak 2 k y K/(l + k 2 ) 2 . This also gives a = for the zonal mode. 

The MHW model spans two limits with respect to the adiabaticity parameter a. In 
the adiabatic limit a —>■ oo (collisionless plasma), the non- zonal component of electron 
density obeys the Boltzmann relation n = n Q (x) exp(ip), and the equations are reduced to 
the Hasegawa-Mima equation 23] . In the hydrodynamic limit a — > 0, the equations are 
decoupled. The vorticity is determined by the 2D Navier-Stokes equation, and the density 
becomes a passive scalar. The advantage of our choice of a as a free parameter is the 



capability for treating the limits in a unified manner. 

The variables in Eqs. (j4j) and (jSJ) have been normalized by 

x/p s ->x, u ci t^t, etp/T e ^<p, ni/n ^n, 

where p s = ^jT c jmuj^ is the ion sound Larmor radius (u s ; = ^/T c /m is the ion sound 
velocity in the cold ion limit), n\ is the fluctuating part of the density. 



In the adiabatic, ideal limit (a = oo, D = 0) the MHW system has two dynamical 
invariants, the energy E and the potential enstrophy W, 

E = \J{n 2 + |Vyf )dz, W=\J{n- Q 2 dx, (6) 

where da; = dxdy, which constrain the fluid motion. Conservation laws are given by 

^ = T„ - D a - D E , ^ = T n - D w , (7) 
dt dt 

where fluxes and dissipations are given by 

f ~ d(p , 
T„ = —k I n—— dec, 



dy 

D a = a I (n — ip) 2 dx, 



D E = D J((V 2 n) 2 + \V(\ 2 )dx, 
D W = D [ (V 2 n - V 2 C) 2 da;. 



Unlike the Hasegawa-Mima model which is an energy- conserving system, the MHW model 
has an energy source T n . Due to the parallel resistivity, n and (p can fluctuate out of phase 
which produces non-zero T n . The system can absorb free energy contained in the background 
density profile through the resistive drift wave instability. 

Note that the same conservation laws hold for the original, unmodified original HW 
(OHW) model, Eqs. (pQ) and (j2J), except that D a is defined by both zonal and non- zonal 
components; _D^ HW = a J (n — <^) 2 dai. In the OHW model, the zonal modes as well as the 
non-zonal modes suffer resistive dissipation. 

We present the linear stability analysis for the zero background (the primary instability). 
Beyond this stability threshold we expect excitation of drift waves. Since the zonal modes 
have linearly decaying solutions, we only consider the form expi(k x x + k y y — uot) (k y ^ 0). 
Linearization of Eqs. (j3J) and (jSJ) around the zero equilibrium (<p = n = 0) yields the 
dispersion relation, 

u 2 + iu(b + 2Dk 4 ) - ibu* - aDk 2 {k 2 + 1) - D 2 k 8 = 0, (8) 

where we defined k 2 = k 2 + k 2 , b = a(l + k 2 )/k 2 , and the drift frequency uo* = k y K,/(l + k 2 ). 
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Solutions to the dispersion relation (jSj) are given by 



U) v = ± 



16^ 

6 2 



cos 



e 

2' 



6 + 2 .DA; 4 =f 6 1 



16a;; 
6 2 



sin 



c<j = dj r + kt?i, tan# = — 4o;*/&. In the limit where Z) = 0, it is readily proved that one of the 
growth rates uj\ is positive if bu* is finite, thus unstable. However, there exists a range of D 
where the drift wave instability is suppressed. The stability threshold is given by 



b + 2Dk A > b 1 + 







sin 



(9) 



6 2 y 2 

and is depicted in Fig. [2j The first unstable mode shown in the figure is the (k x p a , k y p s ) = 
(0,0.15) mode. Below this threshold, an initial perturbation damps out and nothing hap- 
pens. If we choose the parameters in the region beyond the threshold, more than one mode 
starts to grow linearly until the nonlinear terms set in. The left panel shows how many 
modes are excited for given parameters. Most unstable modes are on k x = axis. 
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FIG. 2: Primary stability boundary in cu-k plane and k x -k y plane. 



III. SIMULATION RESULTS 

The HW equations are solved in a doubly periodic square slab domain with box size L = 
2n/Ak where the lowest wavenumber Ak = 0.15 (L ~ 42). The equations are discretized on 
256 x 256 grid points by the finite difference method. Arakawa's method is used for evaluation 
of the Poisson bracket [24j. The time stepping algorithm is the third order explicit linear 
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multistep method 



251 ] . We examine the effects of the parameters k and a on the nonlinearly 



saturated state, and fix D = 10~ 4 throughout this paper. 

We start simulations by imposing small amplitude random perturbations. The pertur- 
bations grow linearly in the initial phase and generate drift waves, then the drift waves 
undergo secondary instabilities which excite zonal flows until nonlinear saturation occurs. 
In the saturated state, we observe that T n ~ D a ^> De,Dw- We compare the MHW and 
the OHW models by showing the spatial behavior of the saturated electrostatic potential 
in Fig. [3j and the time evolution of the total kinetic energy, the zonal component of the 
kinetic energy, and the cross-field transport T n in Fig. HI From Fig. [3] we see that zonally 
elongated structures of the electrostatic potential are generated in the MHW model, while 
rather isotropic vortices are generated in the OHW model. From Fig. H] we see that growth 
of the drift waves is not changed by the modification, but that in the MHW model the 
zonal flows saturate at a higher amplitude (because the modification removes the unphys- 
ical resistive dissipation of the zonal modes). In fact, in the MHW model, the zonal flows 
carry nearly all the kinetic energy in the final state — they have absorbed nearly all the 
energy from the drift waves. In both models, the cross-field transport initially increases as 
the turbulent kinetic energy level increases, but in the MHW model it begins to fall as zonal 
flows absorb the drift wave energy. The build-up of the zonal flow in the MHW model and 
the resulting transport suppression highlight the importance of the difference between the 
MHW and the original HW model in the nonlinear regime 26]. 

Let us show how the parameters k and a affect the saturated state in the MHW model. In 
Fig. El we plot the ratio of the kinetic energy of the zonal flow (F = 1/2 J (d(<p) /dxfdx) to 
the total kinetic energy (E k = 1/2 J \\7(p\ 2 dx) against k and a. It is clearly seen that there 
are two types of saturated states. One is a zonal-flow- dominated state where turbulence 
is almost completely suppressed, and the other is an isotropic turbulence-dominated state. 
The zonal-flow-dominated state suddenly jumps to the turbulent state in a narrow range 
of the parameter space. If we strongly drive the drift wave instability by increasing k, the 
system is likely to reach the turbulent state. From the dependence on a, we can see that 
zonal flows are generated in the adiabatic regime (a ^> 1) while isotropic flows are generated 
in the hydrodynamic regime (a <C 1). These results are compatible with the properties of 
the Hasegawa-Mima model and of hydrodynamic flows as discussed in the next section. 

Let us assume that the generated zonal flows in the y direction can be expressed by a 
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FIG. 3: Contour plot of ip in the saturated state. Zonally elongated structure of the electrostatic 
potential is clearly visible in the modified HW model (a), while isotropic vortices are generated in 
the HW model (b). 



sinusoidal profile, 

V(x) =V sin(Ax). (10) 

The amplitude Vq and wavenumber A = n\Tr/L are determined from the simulation results. 
To estimate A we plot the average wavenumber of the generated zonal flow 

(k x ) = r k( — — . {£ is the kinetic energy spectrum) (11) 

in Fig. [6], and amplitude of the zonal flow in Fig. [71 The average wavenumbers are small 
and rather insensitive to the parameters. This illustrates a feature of 2D flows, which tend 
to generate large scale structures. The wavenumber of a stable zonal flow is typically 0.3 
(corresponding to n\ = 4). The amplitudes of zonal flows are roughly proportional to k, 2 
and are independent of a. 

IV. STABILITY OF ZONAL FLOW 



We examine the stability of the zonal flows obtained from the numerical simulations, 
and compare the stability threshold and the transition point in this section. We consider 
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FIG. 4: Time evolution plots of total kinetic energy, zonal flow kinetic energy and transport of 
MHW and HW models 




FIG. 5: Parameter dependence of the zonal kinetic energy normalized by the total kinetic energy. 
Transitions from a zonal-flow-dominated state to a turbulence-dominated state occur. 



the perturbation around the zonal flow background. The electrostatic potential and the 
density are decomposed as up = ipo(x) + 0(x) exp i(k y y — tot), and n = h(x) exp i(k y y — out) 
where dpo/dx = V gives the background flow in the y direction. By linearizing the MHW 
equations, we obtain an eigenvalue equation containing the effect of k and a, 



dx 2 



kyV 



io 



k y V + ia 



LO 



= 0. 



(12) 
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FIG. 6: Average zonal flow wavenumber versus k and a. 




FIG. 7: Zonal flow amplitude versus k and a. 



We neglect the viscosity. The density fluctuation is determined by 

ia + k v n 



n 



lu — k y V + ia 



(13) 



We solve the eigenvalue equation by the standard shooting method in the domain T> 
{x\ — L/2 < x < L/2}. The boundary is assumed to be rigid <p(±L/2) = for simplicity. 
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A. Hydrodynamic and adiabatic limit 



Before going to the analysis of the HW case, we briefly review the results in two limits: 
the hydrodynamic limit (a — > 0) and the adiabatic limit (a — ► oo). 

In the a — > limit, we recover the Rayleigh eigenvalue equation for neutral fluids, 

~d 2 



dx 2 



k V" 

V UJ — kyV 



= 0. 



(14) 



The well-known Rayleigh 's inflection point theorem demands existence of an inflection point 



for the instability 
Tollmien 



271 ] . The necessary and sufficient condition is also known for this case. 
28| showed the existence of a marginally stable eigenfunction ip s satisfying uj s / k 8j o = 



V(x s ) where x s is the inflection point. ip s satisfies, 



V'l + (A 2 - 



0. 



The solution is given by 
and the critical wave number is 

^s,0 = 



sm{—x) [n : evenj 
cos(^a;) (n:odd) 



(15) 



(16) 



(n = ±l,±2,--0. 



(17) 



If A > tt/L, the marginally stable wave number fc s exists. It should be noted that Tollmien 
does not exclude the possibility that the marginally stable mode is isolated. However, 



perturbation analysis around the marginal 



connected to the marginal solution 29|, I30I ]. 



y mode shows the existence of solutions smoothly 



A similar analysis can be applied to the adiabatic limit, 



dx 2 



UJ — kyV 



= 



(18) 



if k = 0. The marginally stable eigenfunction satisfies, 



+ (A 2 - fc 2 - l)<p B 



0. 



(19) 



The solution is identical with the previous case, but the critical wave number is slightly 
modified to 

»-(£)'-!. (20) 
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FIG. 8: Growth rates for A = 0.3 flow in hydrodynamic limit as described in the text. 

The necessary and sufficient condition of the flow shear for instability is A 2 > (tt/L) 2 + 1. 

We can judge the stability by finding the critical wavenumber. We consider the flow given 
by Eq. (1101) with A = 0.3. The critical wavenumber exists only in the hydrodynamic limit 
for the given profile. On the other hand, the given flow is stable in the adiabatic limit. The 
difference of the two stability conditions (fTT|) and fl20l) comes not from k but from the strong 
coupling between ip and n, and reflects the stabilizing effect of adiabatic parallel electron 
motion. 

Figure M shows the imaginary parts of the eigenvalues for n\ = 4 case in the hydrodynamic 
limit. The eigenvalues are pure imaginary in this limit because of antisymmetry of the flow 
[V(x) = —V(—x)]. Another property in this limit is the scale invariance. The eigenvalues 
do not depend on L and A, but are determined by n\. 

We set Vq = 1. Or, in other words, Vq is normalized out by considering uj/Vq — > to. The 
eigenvalue problem of the given flow profile with ri\ = 4 has the same eigenvalues as that of 
the flow with n\ = 2 in the half domain (solid line). The critical wavenumber for this curve 
is given by fc s ,o( ra A = 2) ~ 0.26. In addition, we find another branch of solutions (broken 
line) which continue to exist until k y < k St0 (n\ = 4) ~ 0.29. 

Next, let us consider the effect of k. Since the critical wavenumber does not exist for the 
profile with A = 0.3, we examine a profile having stronger flow shear by setting L = 5, and 
take n\ = 2 for simplicity. In this setting the marginal wavenumber exists (fc S)0 o ~ 2.14). 
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FIG. 9: Growth rate in adiabatic limit (L = 5, n\ = 2). (a) /c y dependence, (b) k dependence. 

Figure. M shows the eigenvalues obtained in the adiabatic limit for L = 5 and A = 2n/L. 
K is also normalized by k/Vq — * independent of k. Thus the same stability 

condition still holds for finite, but not too large, k. As we see from the figure, the growth 
rate uj\ decreases with increasing k and disappears for large k even though k B>00 exists. We 
need another condition for k. Multiplying ffl8|) by complex conjugate of (p and integrating 
over the domain, we obtain 

r b (v» -l ^ 

(21) 



A I ¥^Mix = 0. 



v \u 



If oj\ 0, V" + k = must be satisfied somewhere in the domain 311] • Applying this 
condition to our assumed flow profile, we obtain the condition k < A 2 for the instability. 
This gives only a necessary condition for the instability, but provides a good estimate [Fig. [9] 
(b)]- 

If we find the eigenvalue u and the corresponding eigenfunction ip, the complex conjugate 
of uj is also an eigenvalue and the corresponding eigenfunction is given by the complex 
conjugate of (p. Thus, we can always restrict our quest for eigenvalues in the upper half 
plane of the complex u plane without loss of generality. This greatly simplifies the situation 
because we can neglect the continuous spectrum on the real u axis. 



B. Hasegawa Wakatani case (intermediate value of a) 



Unlike the previous cases, the complex conjugate of an eigenvalue is not a eigenvalue if 
we include finite a. In this case we must solve for negative cjj as well. Moreover, there exist 
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FIG. 10: Growth rates for HW CclSG clS described in the text. 

two continuous spectra in this case: 

lu = k y V, k y V - ia where \V\ < V . (22) 

Both represent convective transport due to the background flow. One of them is damped 
by the resistivity. These continua may interact with the point spectrum. Thus the situation 
is much more complicated in the intermediate a case compared with the adiabatic and 
hydrodynamic limits. 



First, we show the effect of a and neglect effect of k. We consider n\ = 2 for simplicity. 
Figure \IU\ shows the imaginary parts of the eigenvalues. Three different a cases, and the 
a dependence of the positive branches, are shown. The continuous spectra are shown by 
thick solid lines. In the a = 0.0001 case, two branches from the a — > case (dotted line) 
are also shown for reference, so that it is seen that uj\ is slightly shifted downwards for finite 
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Bifurcation Diagram (D=10 4 ) 
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FIG. 11: Bifurcation diagram showing the correlation between the linearized stability estimates 
described in the text and the regimes observed in our turbulence simulations. 

a. As k y L decreases, the upper (unstable) branch intersects the continuous spectrum at 
marginal stability, and there exists a gap (interval in k y L) occupied by the two continuous 
spectra before this branch continues as a stable mode. The eigenf unctions belonging to the 
eigenvalues in the point spectrum close to this gap become singular. 

For increasing a, we observe the positive eigenvalues disappear at a ~ 0.000417. In 
addition to the two stable branches seen at a = 0.0001, at a = 0.001 another stable branch 
has appeared in the small k y region. By further increase of a we find that the lower two 
branches merge. Beyond this merging point, finite real parts appear, and the eigenmode 
starts to travel in the y direction. 

Our concern is to determine the stability threshold in the a-K plane. Next, we consider 
the effect of k in addition to a. Since k always appears in the form of not and a is small 
in the vicinity of the threshold, the effect of k is rather minor, k does not significantly 
affect the behavior of the eigenvalues except that k controls the amplitude of flow. As we 
stated earlier, the parameters are normalized by Vo, k/Vq k, a/Vo — > a, in the shooting 
calculation, where Vq is proportional to k 2 . 

Finally, we summarize the shooting calculation by showing the bifurcation diagram in 
a-K plane together with the numerically obtained results. The only excitable mode that can 
be resolved in the numerical simulation is the k y = 0.15 mode, which is the first unstable 
mode of the primary instability (see Sec. [TTJ) . In Fig. [TT], we show the stability threshold 
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of ky = 0.15 mode for the primary instability (resistive drift wave instability) and the 
tertiary instability (KH instability). Each mark in the figure denotes a numerically obtained 
saturated state: A, ■. • represent respectively the zonal- flow-dominated, transitional, and 
turbulence-dominated states. In these states zonal flows contain more than 90%, 20-90%, 
and less than 20% of the total kinetic energy, respectively. The qualitative tendency of the 
thresholds in the bifurcation diagram shows agreement between the numerical simulations 
and the KH analysis, i.e. increasing a (k) is stabilizing (destabilizing). Zonal-flow-dominated 
states are observed in between the primary and the tertiary instability thresholds. The 
emergence of a turbulent state is shifted from the primary threshold to the tertiary threshold 
due to the turbulence suppression effect of the zonal flow, which is analogous to the Dimits 
shift observed in ITG turbulence. 

The reasons for the quantitative discrepancy between the boundary of the zonal and the 
turbulent states may be because of the simplification made in the KH analysis; the simplified 
flow profile, the boundary condition and viscosity may also affect the results. 

V. CONCLUSION 

In summary, we have analyzed bifurcation phenomena in two-dimensional resistive drift 
wave turbulence. First, we have performed numerical simulations of the modified HW model 
to study bifurcation structures in a two-parameter (a-n) space. We have shown that, in the 
MHW model, zonal flows are self-organized and suppress turbulence and turbulent trans- 
port over a range of parameters beyond the linear stability threshold for resistive drift 
waves. By performing a systematic parameter survey, we have found that such zonal- 
flow-dominated states suddenly disappear as a threshold is crossed, being replaced by a 
turbulence-dominated state. 

The threshold of the onset of turbulence has been compared with the linear stability 
threshold of an assumed laminar zonal flow profile. Simple theoretical predictions in limiting 
cases explain the qualitative tendency of the stability of the zonal flow, k determines the 
amplitude of the zonal flows, thus, large k destabilizes the zonal flows. On the other hand, 
the adiabatic response of parallel electrons given by a stabilizes them. Numerical analysis 
of the eigenvalue problem determining the stability of the assumed zonal flow profile in the 
HW model also confirms this trend. The constructed bifurcation diagram in the a-K plane 
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for the HW model confirms the scenario of the onset of turbulence in the drift wave/zonal 
flow system being due to the disruption of zonal flows by KH instability. 

The HW model considered here is a particularly simple model, but includes the essential 
physics of interactions between turbulence and coherent structures. This system exhibits 
many other interesting phenomena, but in this paper we have focused on the effect of the 
linear driving term k and the parallel electron response a (including the resistivity). To do 
so, we set the viscosity very small. In this case the zonal flow survives for a very long time. 
However, when the viscosity comes into play, the zonal flows are damped rapidly, and the 
turbulence grows again until zonal flows can be nonlinearly excited and the cycle repeats. 
Thus the system exhibits predator-prey oscillatory behavior. 
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